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ABSTRACT.  The  space  formed  by  the  ground  and 
ionosphere  is  known  to  act  as  a resonator  for 
extremely  low  frequency  (ELF)  waves . Lightning 
discharges  trigger  this  global  resonance , which  is 
known  as  Schumann  resonances  at  the  frequencies  of 
8,  14,  21  Hz  etc.  Even  though  the  inhomogeneity  (like 
day-night  asymmetry,  local  perturbation  etc.)  is 
important  for  such  subionospheric  ELF  propagation, 
the  previous  analyses  have  been  always  made  by 
some  approximations  because  the  problem  is  too 
complicated  to  be  analyzed  by  exact  full-wave 
analysis.  This  paper  presents  the  first  application  of 
the  conventional  FDTD  method  to  such 
subionospheric  ELF  wave  propagation,  in  which  any 
kinds  of  inhomogeneities  can  be  included  in  the 
analysis,  to  be  compared  with  the  observational 
results.  We  show  the  application  of  FDTD  to  our 
problem  and  present  a few  numerical  computational 
results  to  be  compared  with  those  by  the  pre-existing 
analysis  method. 

1.  Introduction 

Recently  there  has  been  again  a great  interest  in 
subionospheric  ELF  waves  because  of  the  two 
important  findings.  The  first  one  is  the  suggestion  by 
Williams  [1]  that  the  global  warming  which  is  a very 
important  issue  for  human  being,  can  be  effectively 


Fig.  1.  Configuration  of  the  Earth-ionosphere 
cavity. 


monitored  by  the  intensity  of  Schumann  resonances. 
Schumann  resonance  is  a global  resonance 
phenomenon  in  the  Earth-ionosphere  cavity  as  shown 
in  Fig.  1,  which  is  triggered  by  lightning  discharges 
in  the  thunder-active  equatorial  region.  The 
resonance  frequencies  are  8,  14,  21  Hz  etc.  in  the  ELF 
range.  The  second  reason  is  closely  associated  with 
the  finding  of  optical  phenomena  (red  sprites  etc.)  in 
the  mesosphere  and  the  lower  ionosphere,  and  it  is 
found  that  this  mesospheric  optical  phenomenon  is  a 
source  of  strong  ELF  signals  (called  ELF  transients) 
[2].  This  ELF  transient  signal  is  one  of  the  important 
tools  for  the  study  of  those  mesospheric  optical 
phenomena  and  then  the  electro  dynamic  coupling 
between  the  atmosphere  and  ionosphere. 

There  have  been  published  a few  monographs 
dealing  with  the  subionospheric  ELF  propagation 
[3-5].  Schumann  [6]  was  the  first  to  predict  the 
presence  of  resonances  in  the  spherical  earth- 
ionosphere  cavity  and  suggested  a mathematical 
formulation  of  the  propagation  problem  at  ELF.  A 
great  simplification  is  the  presence  of  only  a single 
globally  propagating,  zero-order  TEM  mode  [7], 
while  all  the  higher-order  modes  attenuate  severely  at 
distances  of  several  waveguide  effective  heights. 
Despite  this  simplification,  the  complex  structure  of 
the  lower  ionosphere  imposes  an  intricate  three- 
dimensional  electrodynamical  problem  that  cannot  be 
reduced  to  practical  techniques  without  certain 
additional  simplifying  assumptions.  This  is  the 
reason  why  several  fundamental  observed  properties 
of  Schumann  resonances  cannot  be  well  explained 
[5,8],  and  please  look  at  the  forthcoming  monograph 
on  ELF  [8]  indicating  a lot  of  unsolved  problems  in 
the  ELF  range.  The  factors  of  making  the  problem 
very  complete  are  (1)  radially  (vertically) 
inhomogeneity  of  the  ionosphere,  (2)  day-night 
asymmetry,  (3)  local  ionospheric  perturbations  etc. 

In  this  sense  it  is  highly  required  that  very 
complicated  situations  should  be  solved  in  an  exact 
way  even  by  using  the  numerical  methods.  If  the 
exact  numerical  solutions  are  available  for  any 
complicated  cavity,  it  would  be  possible  for  us  to 
examine  the  validity  of  the  previous  approximate 
solutions. 
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2.  Three-dimensional  FDTD  analysis 

The  above  discussion  has  led  to  the  requirement 
for  a three-dimensional  finite-difference  time  domain 
(FDTD)  application  to  our  ELF  propagation  problem 
in  spherical  coodinates  with  azimuthal  dependence. 
Fortunately,  this  3D  FDTD  code  has  been  developed 
by  Holland  [9],  which  can  be  used  in  our  application. 
Fig.  2 (right  panel)  shows  the  spherical  coordinate 
system  (r,  6,  0)  and  location  of  a unit  cell.  The  left 
panel  of  Fig.  2 illustrates  the  location  of  field- 
evaluation  points  on  a unit  cell.  As  shown  in  Fig.  2, 
we  build  a cell  entitled  I,  J and  K along  each  of  three 
axes  ( r , 0,  and  <j>)  . In  the  r direction  1=1  is  the 
starting  grid  (r-  a (Earth’s  radius)  6.4Mm)  and  this 
ground  is  assumed  to  be  a perfect  conductor.  I = Nr 
is  the  outermost  radial  boundary,  is  taken  to  be  r = a 
+ 150  km  and  to  be  a perfect  conductor.  In  the  6 
direction  J = 1 is  the  North  Pole  (0  = 0°)  and  J = Ne 
is  the  South  Pole  (0  = 180°)  . The  coordinate  <p 
indicates  the  longitude  in  such  a way  that  K = 1 and 
K = are  the  same  point,  but  the  latter  has  encircled 
the  globe.  On  this  condition  the  six  electromagnetic 
field  components  are  expressed  by  Holland  [9]. 

The  Maxwell’s  equations  to  be  solved  are 
expressed  as  follows  [10]. 


Fig.  2.  The  spherical  coordinate  system 
{r,6,(j>)  and  location  of  a unit  cell  (right 
panel).  The  left  panel  illustrates  the  location 
of  field-evaluation  points  on  the  unit  cell. 
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(1) 


where  E is  the  electric  field,  H is  the  magnetic  field, 
£0 » are  the  dielectric  constant  and  permeability 
of  free  space,  a is  the  conductivity  which  is  a 
function  of  space,  and  Jex(  is  the  current  source  of  the 
system  (here  this  is  a lightning  discharge).  Here  we 
write  down  the  field  updating  equations  in  the 
spherical  coordinates  as  follows. 


In  the  region  of  1=1  ~NM,  J=2-Ne.3  and  K=2~N^,  the  field  updating  equation  for  Er  component  is  given  by; 
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The  corresponding  equation  for  Ee  is  as  follows. 

E,(UKr«  +fe„/4,  + ff/2F 
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and  this  equation  pertains  to  I-l-N^,  J-2~Ne.i  and  K=2~N^.  And  then,  the  field  updating  equation  for  E{ 
component  is  as  follows. 
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for  I=2~Nr_i,  J=l~Ne_i  and  K=2~N^  Ar,  A0  and  A<j>  are  the  grid  size  in  the  r,  0,  directions,  and  At  is  time  step. 
The  field  updating  equation  for  magnetic  field  components  (Hr,  He,  H^)  are  given  as  follows. 
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for  I=2~N„  J=l~Ne.,  and  K=1~N5,1. 
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for  I=1~Nm,  J=2~Ne_i  and  K=1~NH. 
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for  I=1~Nm,  J=KNe.1  and  K=KN4,1. 

The  field  components  are  on  a modified  Yee  cell  with  unit  vectors  r,  9 , and  q>,  and  the  corresponding 

indices  I,  J,  K.  Spatial  locations  are  given  in  term  of  positions  Ro(I),  0o(J),  <l>o(K).  Points  lying  midway  between 
these  locations  are  given  by  R(I),  0(J),  <j)(K). 


3.  Modeling  of  conductivity  profile  and 
current  source,  and  some  numerical 
examples 

In  order  to  verify  the  usefulness  of  this  3D 
FDTD  computation,  we  show  some  numerical  results 
for  a uniform  cavity  without  any  day -night 
asymmetry,  because  we  can  compare  our  own  results 
with  the  previous  analytical  computations  [11,12]. 


However,  the  height  profile  of  the  conductivity  is 
definitely  induced  in  our  computations. 

First  we  indicate  the  height  profile  of  the 
atmospheric  conductivity,  which  is  expressed  by, 

<j(z)  = a0exp(z!h)  (8) 

where  r is  the  height  above  the  ground  and  h is  the 
scale  height  ( h = 6km).  <J0  is  assumed  to  be 
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10"14[S/m].  This  is  a typical  profile  [10].  Next  we 
assume  a lightning  discharge  current  in  the  form  of 
two  exponentials  by  Bruce  and  Golde  [13],  but  we 
adopt  a model  with  much  slower  time  variation  than 
the  actual  lightning  in  order  to  validate  the  horizontal 
grid  size  of  250km. 

/(/)  = 70{exp(-70/)-exp(-100/)}  (9) 

where  t is  in  second.  Assuming  this  current  source, 
we  made  the  FDTD  computations  in  which  the  grid 
size  in  r is  2km,  and  the  grid  sizes  in  0 and  0are 
/r/72  (about  250km). 

Fig.  3 illustrates  the  computational  results  on  the 
temporal  evolution  of  horizontal  magnetic  field  for 
four  different  distances  (d  = 5,  10,  15  and  20  Mm, 
indicating  the  distance  in  the  6 direction)  when  the 
current  source  is  located  at  the  equator  (0=  nil,  (j) 
~ 0°).  The  pulse  is  generated  at  t = 0 and  we  can 
understand  the  propagation  of  a pulse  in  the  v 
waveguide;  we  notice  the  direct  wave  and  its 
corresponding  antipodal  wave.  As  you  see  from  this 
figure,  with  the  increase  in  propagation  distance,  the 
higher  frequency  components  will  decay  so  that  we 
have  a spread  in  waveform  for  larger  propagation 
distances.  The  horizontal  magnetic  field  of  the 
antipodal  wave  is  opposite  in  sign  to  that  of  the  direct 
wave,  so  that  the  signal  for  d = 20Mm  is  zero  because 
of  the  complete  cancellation  of  the  direct  and 
antipodal  waves.  An  important  characteristic  from 
this  figure  is  that  we  expect  severe  damping  between 
d=5  Mm  and  lOMm,  but  negligible  damping 
between  d~  lOMm  and  15Mm.  These  numerical 
computations  are  compared  with  the  corresponding 
computation  for  their  model  [11]  (the  result  is  shown 
in  Fig.  4),  in  which  there  exist  some  differences  in  the 
modeling.  The  first  one  is  the  difference  in  the  input 
signal  (source)  because  they  used  a delta  function. 
The  sign  of  the  horizontal  magnetic  field  is  taken  to 


Fig.  4.  Computational  result  of  the  temporal 
evolution  of  horizontal  magnetic  field  for 
four  different  propagation  distances  (d=5, 
10,  15  and  20  Mm). 


Fig.  3.  Same  as  Fig.  3,  but  with  the  previous 
analytical  method[ll].  This  figure  should 
be  compared  with  Fig.  3. 

be  oppose  to  that  in  this  paper.  We  have  found  an 
extremely  good  similarity  between  these  two  figures, 
suggesting  the  usefulness  of  this  3D  FDTD  analysis 
in  subionospheric  ELF  propagation.  Another 
difference  is  definitely  due  to  the  difference  in 
propagation  constants.  They  [11]  have  used  a 
heuristic  propagation  constant  (v3  and  v2)  based  on 
the  experimental  data, 

v(/)  = v,-;v2  =(/-2)/6-///100,  (10) 

where  / is  frequency,  and  v3,  v2  are  real  and 
imaginary  part  of  propagation  constant,  respectively. 

On  the  other  hand,  we  have  used  the  atmospheric 
conductivity  profile,  Eq.(8)  as  mentioned  before.  By 
using  the  approximate  formulas  by  Greifinger  and 
Greifinger  [14],  we  can  estimate  the  real  and 
imaginary  parts  of  the  propagation  constant  for  our 
ionospheric  model.  These  results  are  shown  in  Fig.  5. 
Fig.  5(a)  refers  to  the  real  part,  while  Fig.  5(b)  refers 
to  the  imaginary  part  (or  attenuation  constant)  of  the 
propagation  constant.  The  real  parts  are  nearly 
identical  for  the  two  models,  but  we  see  a significant 
difference  in  the  imaginary  part  (i.e.  damping  effect). 
Especially  we  notice  a lot  of  difference  in  the 
attenuation  constant  at  higher  frequencies.  In  the 
paper  [11],  they  used  the  heuristic  propagation 
constant,  Eq.  (10),  so  that  they  have  not  indicated  any 
explicit  conductivity  height  profile.  Of  course  it  is 
theoretically  possible  for  us  to  obtain  the  conductivity 
profile  by  some  approximations  from  Eq.  (10),  but 
this  kind  of  job  is  not  so  meaningful.  In  our  next 
paper  we  will  compare  the  attenuation  from  both 
methods  under  the  same  conductivity  profile  in  order 
to  prove  the  deviation  in  Fig.  5(b)  is  caused  by 
different  computational  models  instead  of  numerical 
errors. 

Finally,  we  show  the  last  figure  for  exact 
comparison  between  ours  and  the  previous  model 
[11]  for  the  same  propagation  distance  (d=10Mm)  in 
Fig.  6. 


Fig.  5.  Comparison  of  propagation  constant,  (a)  refers  to  the  propagation  constant  (real  part),  while  (b)  the 
imaginary  part,  indicating  attenuation  constant.  Two  models  are  compared;  our  model  and 


Nickolaenko  and  Hayakawa  [11] 

4.  Conclusions  and  future  application  of 
this  FDTD  analysis  to  much  more 
complicated  models 

In  the  previous  section  we  have  shown  only  one 
example  of  our  FDTD  analysis  for  subionospheric 
ELF  propagation  and  a comparison  with  the  previous 
analytical  solution  has  indicated  that  our  FDTD 
application  would  be  of  great  potential  in  our  future 
studies.  As  you  see  from  the  above  discussion,  the 
field  updating  mechanism  is  exactly  the  same  as  that 
of  the  high-frequency  FDTD  [16],  and  it  seems  that 
there  is  no  special  treatment  in  the  ELF  wave  analysis. 
We  can  list  what  we  will  be  able  to  do  in  this  ELF 
field  by  means  of  our  FDTD  analysis. 

As  for  the  ELF  propagation  itself  (such  as  the 
propagation  of  ELF  transients  (Q  bursts,  slow  tails)), 
we  can  study  the  reflection  mechanism  of  these  ELF 
waves  in  the  lower  ionosphere  in  details  by 
comparing  our  FDTD  results  with  the  previously 
proposed  approximations  [14,15]. 

Then,  as  for  Schumann  resonance,  there  are 


Time  (ms) 

Fig.  6.  Comparison  of  the  ELF  waveforms  for  a 
particular  d=10Mm  estimated  by  our 
method  and  previous  method  [11]. 


several  possibilities  of  the  application  of  our  FDTD 
method.  For  example,  (1)  die  full  consideration  of 
day -night  asymmetry  in  order  to  well  explain  the 
diurnal  variation  of  Schumann  resonance  intensity 
and  frequencies,  (2)  the  local  perturbation  (such  as 
the  perturbation  at  high  latitudes  (auroral  regions)). 

We  will  use  our  FDTD  analysis  for  these  above 
questions  in  order  to  improve  the  physics  coming 
from  these  subionospheric  ELF  propagation  studies. 
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